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Abstract 


This paper presents a modification of the spring analogy scheme which uses axial linear spring stiffness with selective spring 
stiffening/relaxation. An alternate approach to solving the geometric conservation law is taken which eliminates the need for storage 
of metric Jacobians at previous time steps. Efficiency and verification are illustrated with several unsteady 2-D airfoil Euler 
computations. The method is next applied to the computation of the turbulent flow about a 2-D airfoil and wing with two and three- 
dimensional moving spoiler surfaces, and the results compared with Benchmark Active Controls Technology (BACT) experimental 
data. The aeroelastic response at low dynamic pressure of an airfoil to a single large scale oscillation of a spoiler surface is computed. 
This study confirms that it is possible to achieve accurate solutions with a very large time step for aeroelastic problems using the fluid 
solver and aeroelastic integrator as discussed in this paper. 


1. Introduction 

Aeroelasticity is the study of the interaction of inertial, elastic 
and aerodynamic forces acting on a structure. 
Aeroservoelasticity includes the additional mutual interaction 
of a control system and the aeroelastic structural response. 
When constructing computational algorithms to model 
problems in these disciplines, accuracy and robustness must 
receive due attention. This is especially true in view of 
recent interest in the accurate computation of fluid-structure 
interaction in the presence of a strongly nonlinear flow field. 
This includes examples such as the computation of limit 
cycle oscillation and the aeroelastic response of a structure 
due to large-scale control surface motion. To address these 
issues of robustness and accuracy, recent aeroelastic research 
has focused on the development of new algorithms that 
integrate the structural equations of motion in 
synchronization with the aerodynamic equations of motion. 
Among the methods used or promoted are the closely 
coupled lagged approach 1 ' 2 to updating the structure 
equations and the implicit iterative coupling of both structure 
and fluid. 5 Other approaches recently discussed include the 
Arbitrary Lagrangian-Eulerian (ALE) J " < and the Implicit 
Continuous-fluid Eulerian (ICED-ALE) methods. 6 Recent 
mesh deformation algorithm development has emphasized 
the robustness and efficiency necessary for coupled structure- 
fluid time marching computations. 7 
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The purpose of the present paper is to revisit old approaches 
to mesh deformation and the integration of the structural 
equations of motion. In particular the spring analogy mesh 
scheme is revisited with modifications that enhance its value 
for structured grids with complex geometry. The finite 
dimensional state-space predictor-corrector method 8 also 
merits renewed consideration in view of recent advances in 
computational fluid solvers. The performance of the 
modified spring analogy and the predictor-corrector 
integrator will be assessed in this paper. 

In the development of the present method, several 
modifications of the typical approach to CFD with a 
deforming mesh have been made. For unsteady problems, a 
self-consistent approach to the Jacobian of the coordinate 
transformation is to compute it by solving an independent 
equation called the geometric conservation law (GCL). 59 ' 10 
In reference 3, the GCL fluxes are retained in the Navier- 
Stokes equations as source terms, replacing the time rate of 
change of the volume. Not only does this result in self- 
consistency, it also eliminates the need for storage of either 
the discrete time derivative of the Jacobian or of the Jacobian 
itself at several previous time steps. A variation of that 
approach is used here. The metric fluxes are retained in the 
present scheme on both the left and right hand sides of the 
discrete approximate factorized (AF) equations in a cell 
centered finite volume discretization of the equations. In 
reference 3, the flow equations were solved in finite 
difference form using central differencing of the flux terms 
with explicit and implicit artificial dissipation. Several 
options are also available in the present code for construction 
of fluxes and time accuracy. Third order upwind Roe's flux 
differencing and second order backward time differencing are 
used here. Subiterations either include the physical time step 



(t-TS subiteration) only or both physical and “pseudo” time 
step based on local Courant-Friedrichs-Lewy (CFL) number 
for convergence acceleration (t-TS subiteration). This affords 
the opportunity to assess the effect of the increased 
robustness of a higher order upwind scheme with pseudo 
time subiterations on the simulation of unsteady aeroelastic 
response to a nonlinear flow field. Finally the present 
deforming mesh scheme is used with a time accurate field 
equation turbulence model. 

Another aspect of this work is the modification of the spring 
analogy scheme with axial stiffness for use with a structured 
grid. References 6 and 11-12 give recent examples of its use 
with structured and unstructured grids. There are advantages 
to the spring analogy, most notably its simplicity and ability 
to smooth grids. However, it does not ensure positivity of 
control volumes, and in its most typical form is memory 
intensive. The modified spring analogy using nonlinear 
torsion stiffness has been shown to preclude non positive 
control volume (and thus the more restrictive case of grid 
point crossing). 7 It does not alone allow spacing control and 
does not include the smoothing available from the axial 
spring analogy. Thus it would seem that an optimal approach 
in the context of the spring analogy would be a combination 
of several stiffness types. It will be shown that in its present 
form the spring analogy scheme with linear axial stiffness is 
efficient and with a structured grid also allows some savings 
in memory. It can therefore provide a starting point for a 
more elaborate deforming mesh scheme. It is not, however, 
without problems. Non positive grid volumes can still be 
encountered. A partial solution is presented that is useful for 
at least mildly complex geometry, but it must be emphasized 
that the present work is considered the initial step in the 
development of a robust deforming mesh scheme. The 
intention is that future developments will include the 
capability of grid adaptation and torsion stiffness that will 
allow more complex geometry and control over grid 
orthogonality. 


capability to dynamically oscillate a grid, 18 a feature that will 
be used in the code comparisons that follow. Since the form 
of the thin layer Navier-Stokes equations as used here and the 
solution procedures are not typical, the equations and 
details of the solution will be presented. In a generalized 
deforming coordinate system the differential form of the 
equations can be written 
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and Q = (p,pu,pv,pw,e) ■ The inviscid flux and 
viscous diffusive terms are (f.G.h) and ( F r ,G v ,H v ) 
respectively. Note that the time derivative on the left-hand 
side of equation 2.1 is not written in strong conservation-law 
form, and that the last bracketed item in the residual is 
composed of the grid speed fluxes arising from the GCL. 
Letting 
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The scheme is first verified for several unsteady 2-D airfoil 
Euler computations and several unsteady oscillating spoiler 
Navier-Stokes computations using the Spalart-AUmaras 
turbulence model. The oscillating airfoil and spoiler 
computations explore time accuracy and geometric modeling 
issues. Among other features, the oscillating spoilers are 
modeled using rigid body rotation rather than the more 
common shearing control surface motion. This adds 
somewhat to the realism of the model. A steady rectangular 
wing 3-D spoiler case and an unsteady 3-D oscillating spoiler 
case have been computed and are compared with 
experimental data from the Benchmark Active Controls 
Technology (BACT) test. 13 ' 15 The BACT test data presents 
steady, unsteady and aeroelastic pressure and flutter onset 
data from a test of a rectangular NACA 0012 wing conducted 
in the NASA Langley Transonic Dynamics Tunnel. Data 
associated with that test is available for an active flutter 
suppression control system using spoiler and trailing edge 
control surfaces. This present work is the initial part of an 
effort to utilize this data for the validation of a CFD, 
aeroelastic, control system simulation methodology that 
allows for moderately complex control surface motion. 

11. Thin layer Navier-Stokes equations 

The finite volume thin layer Navier-Stokes code that forms 
the basis for this new development is CFL3D version 5.0. 
This code has been used extensively to simulate many 
unsteady 2 and 3-D single and multizonal problems. 16 ' 17 It 
includes multigrid and grid sequencing capability and has 
available many choices for turbulence model. It also has the 


and constructing the remaining inviscid and diffusive 
fluxes, the complete equation set is obtained. In integral 
form, the equations 2. 1 can be written 

l \ v J-'Q,dVdT = (E-gQ]dAdr 

- jJ v (V E v )dVdT 
-[\ v Q(J-'Vg)dVdr (2.3) 

where g is the grid speed, E = (F,G,H) and 

£ = , G , H ) , while f 1 is the ratio of the physical to 

computational volumes. The fluxes are evaluated at the cell 
face centers, while volume integrals over the computational 
cell volumes V are evaluated discretely at control volume 
centers, making this a cell centered finite volume scheme. 

The residual composed of the right hand side of equation 2.3 
is linearized about the most recent subiteration m and the 
resulting equations are approximately factored. The resulting 
discrete equations take on the form 
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biased Roe’s flux difference splitting and a min-mod flux 
limiter. Multigrid is used in the following computations. 
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for the residual as defined above and 
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In the equations above M = dQ / dq , 
A’ = 3 (F — F v )/dq and so forth. The solution variable is 


q=(p,u.v.w,p). Since the fluxes contained in equation 2.5 
and 2.2 (last term) are similar in form to the fluxes of the 
flow quantities both can easily be computed within the same 
subroutines. This has required the addition of 20-40 coding 
lines to each of the inviscid flux subroutines (right and left- 
hand side terms). Furthermore, the desired first or second 
order accuracy in time is automatically achieved by making 
use of the metric time derivatives already being used in the 
spatial fluxes of conserved flow quantities. Since the 
volumes at several previous time steps are not stored, very 
little additional memory and no change in the overall code 
data structure is required. Using this approach does require 
that equation 2.5 and the metric fluxes of equation 2.2 be 
recomputed at each subiteration and added to the equations 
2.4a-c, resulting in what appears to be a modest added 
computational effort. 


The metric fluxes in the new terms are evaluated at cell face 
centers, and are also differenced in time in a way that is 
consistent with the temporal differencing of conserved flow 
variables. This results in consistency between new metric 
fluxes of equation 2.5, the metric time derivatives in the 
physical fluxes and the integral of J 'Q, over the control 
volume, and finally in a temporally and spatially consistent 
solution of the GCL. The value of the metric Jacobian at 
time step n + / is used everywhere in equations 2.4 and 2.5. 
This and the evaluation of the residual at iterate m implies a 
time discretization of the volume integral equation 2.3 
centered about time steps n + 3/2 and n + 1/2. The equations 
2.4 also include local CFL based "pseudo" time stepping (I- 
TS subiteration) for convergence. Depending on the values 
of 0 and <jJ , the option of first or second order physical and 
pseudo time stepping can be exercised. In the present 
computations the temporal discretization of equation 2.3 is 
accomplished via second order backward differencing (i.e. 0 
= 1/2). The discrete inviscid fluxes use third order upwind- 


CFL3D version 5.0 includes many turbulence models. The 
turbulence model used in the present computations is the 
Spalart-Allmaras field equation model, selected here 
primarily because of its robustness at large time steps. Since 
the differential equation for eddy viscosity in the Spalart- 
Allmaras model is solved in nonconservation form it does not 
impose any additional requirements on the metrics due to 
grid deformation. This approach simplifies the procedure 
when using unsteady field equation turbulence models, since 
for unsteady computations no additional treatment of metrics 
or alteration, other than consistent time differencing, is 
required by the solution of the turbulence model. 


111. Mesh scheme 


The mesh scheme is a modification of the spring analogy 
using axial spring stiffness. The approach of Batina can be 
written 


MS+CS+KS= 0 (3.1) 

where K is the spring stiffness matrix, M and C are the mass 
and damping matrices and the mesh displacement is given 
by 5. In this formulation M = C = 0. The stiffness would be 

-II 1/2 . . 

r J- r lh and ' 

and j represent two adjacent nodes. For an unstructured grid, 
storage of stiffness values based on the current locations of 
each node pair is required. The spring stiffness would be 
updated as the mesh deforms. For a structured grid the 
problem is somewhat simplified. Spring stiffness in the mesh 
interior can be controlled by the spacing of the appropriate 
boundary grid points. For instance, if the nodes at (i,j,k) and 
( i,j,k+I ) are considered and if boundary spacing at 
boundaries i = 1 and i = i, mx are used 
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where m=k+l designates in this case the volume edge 
between the k and k+I grid points, and 
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Since these stiffness values are set at the start of the 
computations they do not vary in time. They also require 
storage only for each of the six computational boundaries. 
Linear spring stiffness allows the possibility of grid point 
crossings, but in view of the fact that this is a restricted case 
of the general problem of grid lines crossing, the better 
solution is the addition of nonlinear torsion stiffness. 7 This 
can be added to the present spring stiffness at a future time. 
Finally, it is commented that the axial stiffness approach used 
here results in smoothing of the mesh and also allows 
adaptation based on the flow solution. 



The problem of grid collapse around convex surfaces is 
handled by selectively increasing/decreasing stiffness based 
on surface curvature. Stiffness values in two coordinate 
planes normal to a surface are varied based on surface 
curvature in the coordinate projection of those planes onto 
the surface. The final mesh is the weighted combination of 
the two planar solutions. Take for example a ^r| -plane solid 
surface at k-l from which curvature information is required 
to construct the grid. Mathematically the expression for 
interior grid point location can be written 

n jk =(1“ £)[/, X + /: X 
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In equation 3.2 k ml —\k m e C|A ’” and 

k m -, — -7 k m e C> ' 1 for the volume edges including the 

endpoints k+l and k-l, and k ml — k m -, = k m otherwise. 
The index ml in the summation ranges over the grid points 
adjacent to ijk in the 4C-plane (indices j.k) and index m2 in 
the qij-plane (indices i,k). The constant C/ is a gain factor 
that adjusts sensitivity to surface curvature. A value of 
Ci-20 seems to work well for many geometries. The surface 
curvature parameters A5 and A n can be arrived at by 
considering, for example, the surface grids of Figure 1 in the 
\ O' index) direction. A measure that accounts, at grid i.j.k, 
for convex or concave curvature at the surface point i.j and 
k=l is 
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In terms of grid geometry this is 
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defines an outward normal to the surface. Setting 
£ = exp(-C 2 &) ensures that at least several grid points 
near the surface remain orthogonal to the surface. 

IV. Aeroelastic method 

The time marching simulation of the aeroelastic responses is 
obtained using the finite dimensional state space predictor- 
corrector method, described in references 8 and 19-20, that 
solves the de-coupled modal structural equations of motion. 
Several versions of the finite dimensional state space variable 
method are presented in those references. Of those the 
predictor-corrector method offers the best performance. The 
predictor step marches the structure using the solution of the 
modal equations at step n to get the surface deflection at time 
step n+l. This solution is based on a second order 
(trapezoidal) extrapolation of the generalized aerodynamic 
forces at n and n-l. This provides the surface shape for a 
recomputation of the fluid mesh and the fluid domain 
solution at n+I. After the solution of the fluid domain, the 
corrector step then solves the modal equations at the time 
step n+I using the averaged generalized forces at n and n+l. 

V. Results 


The first computational results illustrate robustness and 
efficiency. A comparative study was made for a pitching 
airfoil using the present deforming mesh scheme in which 
only the inner boundary is oscillated and that in which the 
entire grid was dynamically oscillated. The capability to 
dynamically oscillate the (rigid) grid has been developed and 
verified in a separate paper. 18 Both mesh movement schemes 
used the same grid having dimensions 297x105. Test cases 
were the AGARD NACA 0012 cases 3 and 5. 21 In each, the 
number of subiterations and the t-TS local CFL number 
were adjusted to give a convergence of the Lj-norm of at 
least 10' 7 at each time step. The results are virtually 
indistinguishable, as seen in Figure 2 for the Euler 
computation of AGARD Case 5. This is true even at 16 t.s. 
per cycle, for which At=3.2 (nondimensionalized by speed of 
sound), requiring between 20-30 subiterations and some 
tuning of the local CFL number to adequately converge. An 
additional test of the time step convergence using the 
deforming mesh capability is also shown in Figure 2. The lift 
coefficient shows no discemable variation, while the moment 
coefficient, although overall nicely convergent, reveals an 
anomalous lack of convergence around zero degrees (upward 
motion). At 128 time steps per cycle, any difference in the 
deforming mesh and rotating frame solutions is imperceptible 
with the exception of moment coefficient on the downward 
side of its motion at around -2 degrees. (Figure 3) In Figure 
4 a comparison of Navier-Stokes solutions at 64 and 2048 
time steps per cycle of airfoil oscillation for AGARD case 3 
is presented. This case involves slight trailing edge 
separation. This, and all the computations to follow, are 
made with the Spalart-Allmaras turbulence model. As can be 
seen in the figure, very slight differences between the 64 
t.sVcycle solution and that at 2048 t.s./cycle only show up 
during trailing edge separation just after the start of pitch 
down. 


In practice the values of At and A n are limited in upper and 
lower bound, currently At=min(C|At,0)/Ci or A ; =max(C|A t ,- 
12)/Ci and similarly for A n . An update of/, (A$) and f 2 ( A„) 
are required once per time step for each surface or grid edge 
plane. This leaves only the computation of equation 3.2 to 
solve for the interior of the mesh. This is a linear equation 
typically requiring 1-2 iterations of a Gauss-Seidel procedure 
for moderate motions of the grid. Finally, the function P 


The preceding Euler computations afford the opportunity to 
assess performance. The grid deformation scheme, the added 
terms for volumetric change in the Euler equations and the 
recomputation of cell volumes combined add slightly less 
than 4 percent CPU time (128 t.s./cycle case, 6 subiterations) 
over that of the unsteady code with rotating frame. This is 
with the grid, volumes and metrics updated at each time step. 
The percent of CPU time is. of course, reduced if larger time 



steps and more subiterations are used. Total code 
performance is approximately 1 1.83 p sec/cell/iteration C-90 
time for the two-dimensional Euler computations. 

The final results of this paper are aimed at assessing what is 
required to compute steady and unsteady two and three 
dimensional spoiler flows. Specifically, the issues addressed 
are the geometric modeling of 2-D and 3-D spoilers and the 
convergence and time step behavior of the present fluid and 
aeroelastic solvers. In each of the following cases the 
airfoil/wing-spoiler geometry are based on that used in the 
Benchmark Active Controls Technology (BACT) NACA 
0012 test. 11 ' 15 This allows a comparison of the present results 
with the BACT experimental database of oscillating spoiler 
cases. The spoiler is modeled here as a ramp and backward 
facing step. A two dimensional version is shown in Figure 5. 
The three dimensional version, seen in Figure 7, is modeled 
as a backward facing wedge having finite spanwise width. 
As can be seen in both figures, the backward step has a slight 
slope modeled by spacing three grid points between the 
spoiler trailing edge and the wing surface. This clearly does 
not model the cavity beneath the spoiler, nor the gap between 
the spoiler and flap leading edge. As the spoiler moves the 
actual spacing over the backward step is constrained via arc 
length to remain more or less true to the original grid spacing 
over the wing. Spoiler motion is modeled with rigid body 
rotation about the spoiler hinge line. The typical approach to 
modeling the spoiler deflection is to shear the wing surface 
grids or, at best, to use a two coordinate modal control 
surface motion model, either of which result in surface 
warping. The present approach embodies true control surface 
motion. 

Before comparing computed results with the steady and 
unsteady experimental data, unsteady Navier-Stokes 
solutions (and grids for a 2-D configuration) are presented in 
Figures 5 and 6. This computation investigates time step 
convergence, which is important before launching into the 
more costly 3-D oscillating spoiler computations. These 
computations are especially challenging since they involve 
the growth and collapse of a large region of recirculating 
flow. They require significantly more subiterations (50-60 at 
64 time steps per cycle) and tuning of the x-TS CFL (for 
extremely large physical time steps) to reach reasonable 
accuracy at each time step. Nevertheless, the convergence at 
each time step was carefully monitored. The Lj-norm was 
maintained at approximately 10‘ 7 or better at each time step. 
As can be seen in Figure 6, there is no discemable difference 
in the lift and moment coefficients in going from 1024 to 64 
time steps per cycle of spoiler oscillation. However, the out- 
of-phase part of the transform of the unsteady pressure 
coefficient shows some reduction in accuracy for the larger 
time steps. While all the computations appear to be 
sufficiently accurate for engineering purposes, it is clear that 
for high accuracy, 260-300 t.s./cycle or more are required. 
This is true, say, for simulation of response to control system 
input at a condition near flutter onset, which is typically 
sensitive to the phase angle of the aerodynamic forcing. 

A comparison of steady and unsteady computed results will 
now be made with experimental data. The purpose here is to 
assess the computational model for a surface geometry 
having moderate complexity. The cases computed are for 
both a 3-D steady and oscillating spoiler on the rectangular 
NACA 0012 wing used in the Benchmark Active Controls 
Technology (BACT) Test. As discussed earlier, surface 
geometry and grids are shown in Figure 7. The deflection in 
the steady case is t, e = 15 degrees. In all the other 
computations viscous terms were included only in the 
direction norma! to the solid surface. It was found after the 


fact that the flow around the trailing edge of the spoiler is 
moderately altered when streamwise viscous terms are 
included. For this reason the steady 3-D spoiler calculation 
included viscous terms in all three index directions. 
Resolving the grid at the trailing and spanwise edges of the 
spoiler and the addition of the streamwise viscous terms 
results in a better match with experiment near the spoiler 
trailing edge. 

The grid for the steady case has somewhat finer resolution 
than that used in the unsteady 3-D computations. It has 
dimensions 201x73x73. Steady pressure coefficients at span 
locations of 40 % and 60% are shown in Figure 8. These are 
just inboard and at mid span of the spoiler which has 
spanwise edges at 45% and 75% span. The pressure 
distribution in the region of the spoiler and its wake show 
quite good agreement with experiment. One minor exception 
is at 60% span where a slight divergence is observed just 
ahead of the spoiler trailing edge and just ahead of the shock. 
The solution at the 40% span location is off somewhat from 
.6c to ,75c. This may be due to an inadequate geometric 
modeling of the spoiler side edges resulting in an improper 
side wash and vortex development. Figure 9 compares the 
computed real and imaginary unsteady pressure coefficients 
for an oscillating spoiler case with experiment. The grid used 
for this computation was 153x53x41 in dimension. The 
mean spoiler angle is 5 degrees and the amplitude of motion 
is 4.5 degrees. This was computed with 128 t.s./cycle for a 
reduced frequency based on semi-chord of k=0.1088 (9.56 
Hz). An obvious notable difference between the 
computations and experiment is observed in the real part of 
the unsteady pressure coefficient at 60 % span. While not 
verified, experience with steady cases indicates that some of 
this difference may be due to the thin-layer viscous 
approximation used. The other significant difference is in the 
imaginary part of the unsteady pressure coefficient around 
the trailing edge of the spoiler at 60% span. 

The simulation of an aeroelastic transient induced by a single 
oscillation of a 2-D spoiler for an airfoil subject to pitch and 
plunge is presented in Figures 10-12 and Tables 1 and 2. 
This case provides a nominal test of the performance of the 
complete fluid-mesh-structure solution technique. The 
aeroelastic response is due to a single large-scale asymmetric 
oscillation of the spoiler. The spoiler was oscillated at 
approximately the free vibration frequency of the plunge 
mode. The dynamic pressure = 20 psj) is low enough to 
give a moderate aeroelastic response for this size of spoiler 
motion. After an aeroelastic response using the mesh shown 
in Figure 5 was computed, it was found that the aeroelastic 
transient response was sensitive to grid spacing asymmetry 
between the upper and lower surfaces. A mesh that is 
symmetric and has the same spacing over both the upper and 
lower surface was used as the final grid. This grid had 
dimensions 319x105. Furthermore, the quality of the grid 
especially around the spoiler region was found to 
considerably impact the time step convergence of the 
aeroelastic solutions. After some effort at improving grid 
quality, the described solutions were obtained. 

The solutions of Figure 10 give transient pitch and plunge 
response calculated with four time step sizes. Note that the 
maximum CFL number for the largest time step (At-. 5859) 
is approximately 2xl(f. The convergence of the plunge mode 
response with time step size reduction can be clearly seen in 
the detailed inset figure. Table 1 summarizes data from these 
runs. The fluid domain subiterations shown in this table are 
what were required to converge the solution (L> norm of 
density) to around I0' 7 at each time step. Between 2048 and 
8192 t.s./plunge cycle there is virtually no difference in the 



damping ratios. Between the largest and smallest time steps 
the difference is less than 1%. The convergence of the 
aeroelastic response with number of subiterations at the 
largest time step can be seen in Figure II. It appears from 
this figure that 32 subiterations of the fluid domain per time 
step represent a very nearly converged solution. One can 
also see from this figure and from Table 2 that if less 
accuracy is required, fewer subiterations will also result in a 
solution. On the other hand, a computation that was 
performed using only a physical time step (t-TS subiteration) 
was unstable for time steps larger than At-.OOl. In fact, the 
solution at At=.0012, shown in Figure 12, diverged shortly 
after the final time step shown. Even if it had not diverged, 
this small time step makes this problem solution infeasible 
with the t-TS method, and necessitates CFL based local time 
step subiteration (identified as t-TS in the figure). The 
accuracy of the t-TS solution is clearly not established with 
only one computation. However, if the solutions with t-TS 
subiterations are converging to the correct solution, as 
appears from the figure, the t-TS solution even at At =.00 1 
does not appear to be more accurate than the t-TS solution at 
a time step 100 times larger. 

VI. Concluding remarks 

A numerical scheme for the Euler and Navier-Stokes 
equations in a general dynamically deforming coordinate 
system that incorporates the geometric conservation law into 
the governing equations has been evaluated in this paper. By 
casting the equations in the form shown here, a somewhat 
more economical scheme is obtained for problems requiring 
deforming meshes. A structured mesh scheme, based on the 
spring analogy, has been applied to moderately complex 
surface geometry. Its performance has been assessed for 
time steps spanning a large range, applied to problems from 
oscillating two-dimensional airfoils and spoilers to a wing 
having a three-dimensional oscillating spoiler surface. The 
spoilers have been modeled as a ramp and backward step and 
the oscillation by true rigid body rotation. The mesh scheme 
does not appear to result in mesh entanglement for the 
geometry and mesh motion considered here even with 
extremely large time step {CFL^2xlCf). This is due both to 
the fact that the fluid mesh near the moving surface 
undergoes motion prescribed by the surface and to the 
selective relaxing/stiffening of the spring stiffness in the grid 
deformation scheme. Even at very large time steps the 
scheme appears to retain considerable accuracy if a sufficient 
number of t-TS subiterations are used. In the case of 
strongly nonlinear flow driven by spoiler oscillation, 
numerical results have shown that the error due to increasing 
time step that does arise is mainly in the out-of-phase 
component. This has implications for the accurate 
computation of flutter onset and other aeroelastic 
phenomena. The steady and unsteady computed results for 
the three-dimensional spoiler and wing qualitatively compare 
well with experiment, with the exception of the fairly small 
out of phase component of the unsteady pressure coefficient. 
Finally, aeroelastic capability has been demonstrated for the 
damped transients of a wing having pitch and plunge modes. 
The case computed was at low dynamic pressure (q„ = 20 
psf) in which the modes were initially perturbed by a single 
large-scale oscillation of the spoiler surface. The results show 
that accuracy at very large time step sizes, similar to that 
used in the oscillating spoiler computations, can be achieved 
here as well. This suggests that for these cases and solution 
method, the convergence of the fluid solution is the dominant 
factor in the accuracy of the structure-fluid solution. Again, 
this is possible because of the t-TS flow field solution. 
Because of stability limitations, the largest time step possible 


with t-TS subiteration was at least 600 times smaller than the 

largest time step with t-TS subiteration. 
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Figure 1 . Orientation of surface and interior grids. 
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Figure 2. Time step study for deforming mesh scheme. 
AGARD NACA 0012 Case 5, NACA 0012, M„ =0.755, 
a = 0.016 + 2.51 sin(2ktM„) degrees, k = 0.081 




Figure 3. Comparison of results for rotating frame and deforming mesh. 
AGARD NACA 0012 Case 5, NACA 0012, M„ =0.755, 
a = 0.016 + 2.51 sin(2ktM,J degrees, k = 0.081 




Figure 4. Time step study for deforming mesh viscous computation. 
AGARD NACA 0012 Case 3, M_ = 0.60, a = 4.86 + 2.44 sin(2ktMJ 
degrees, k = 0.081 
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Figure 5. Two-dimensional mesh for oscillating spoiler study. 




Figure 7. BACT model upper surface and mesh. 
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Figure 8. Comparisons of BACT model steady pressure coefficients. 
M,„ =0.77, a = 0 degrees S sp = 15 degrees. 
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Figure 9. Comparisons of BACT model unsteady pressure coefficients. 
= 0.77, a = 0 degrees 8 sp = 5+4.5 sin(2 MJct) degrees. 
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Figure 10. Computed aeroelastic response due to a 
single spoiler oscillation (9.5 degree excursion), 
2D BACT model. M„ = 0.77, a = 0 degrees, 
q„ = 20 psf, U„ = 373 fps. Figure 10). 


Table 1. Computed damp- 
ing ratios for two aeroelastic 
modes, (Same model and 
condition as Figure 10). 





Figure 1 1 . Computed aeroelastic response 
due to a single spoiler oscillation. (Same 
model and condition as Figure 10, 
At=0.5859.) 
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Table 2. Effect of number of 
subiterations on damping ratio, 
(Same model and condition as 
Figure 10) 
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Figure 12. Comparison of t-TS and t-TS time stepping for aeroelastic 
response (same model and condition as Figure 10). 




